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Abstract 

We study a simple model of a self-assembled, room temperature Coulomb-blockade 
nanostructure containing a metallic nanocrystal or grain connected by soft molec- 
ular links to two metallic electrodes. Self-excitation of periodic grain vibrations at 
10 - 100 GHz is shown to be possible for a sufficiently large bias voltage leading to 
a novel 'shuttle mechanism' of discrete charge transfer and a current through the 
nanostructure proportional to the vibration frequency. For the case of weak elec- 
tromechanical coupling an analytical approach is developed which together with 
Monte Carlo simulations shows that the shuttle instability for structures with high 
junction resistances leads to hysteresis in the current - voltage characteristics. 

Key words: Mesoscopic physics, Coulomb blockade, self-assembled structures, 
electron tunneling, micromechanics 



1 Introduction 



Conventional microelectronics is approaching a limit where further miniatur- 
ization is no longer possible. This has motivated a vigorous search for alter- 
native technologies such as 'single electronics', which is based on Coulomb 
charging effects in ultrasmall structures [1,2]. A novel approach to building 
such structures — nature's own approach — is self-assembly using molecular 
recognition processes to form complex functional units. Familiar examples of 
molecules with this ability are amphiphilic chain molecules (e.g. thiol) in solu- 
tion which form ordered films on easily polarizable metal surfaces, antibodies 
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which find and bind to specific molecular targets, and DNA strands which rec- 
ognize and bind to matching sequences. Recent progress include the successful 
use of self-assembled DNA templates for making tiny silver wires connecting 
macroscopic gold electrodes [3] and the demonstration of room-temperature 
Coulomb blockade behavior in novel composite mesoscopic structures contain- 
ing both metallic elements and self-assembled organic matter [4-6]. Charging 
effects in the latter structures are in the focus of the work we present here. 

The crucial aspect of the new room-temperature Coulomb blockade structures 
from the point of view of our work is that they contain metallic grains or molec- 
ular clusters with a typical size of 1-5 nm that can vibrate; their positions are 
not necessarily fixed. This is because the dielectric material surrounding them 
is elastic and consists of mechanically soft organic molecules. These molecular 
inter-links have elastic moduli which are typically two or three orders of mag- 
nitude smaller than those of ordinary solids [7] . Their ohmic resistance is high 
and of order 10 7 - 10 8 ohm, while at the same time they are extremely small 
- a few nanometers in size. A large Coulomb blockade effect in combination 
with the softness of the dielectric medium implies that charge transfer may 
give rise to a significant deformation of these structures as they respond to 
the electric field associated with a bias voltage. Hence the position of a grain 
with respect to, for instance, bulk metallic leads is not necessarily fixed. For 
the model system of Fig. 1 — containing one metallic cluster connected by 
molecular links to two metallic electrodes — we have recently shown [8] that 
self-excitation of mechanical grain vibrations at 10 - 100 GHz accompanied by 
barrier deformations is possible for a sufficiently large bias voltage. This effect 
amounts to a novel 'shuttle mechanism' for electron transport. 

The purpose of the present paper is to carry the analysis presented in Ref. [8] 
quite a bit further. Strictly speaking the analytical part of the analysis in 
our earlier work is valid for low tunnel-barrier resistances when the rate of 
charge redistribution between grain and leads, inversely proportional to the 
tunneling resistance, is assumed to be so large in comparison with the vibra- 
tion frequency that the stochastic fluctuations in grain charge during a single 
vibration period are unimportant. In order to describe the opposite limit of 
low charge redistribution frequencies characteristic of high-resistance tunnel 
barriers we develop here a new approach for considering the coupling betwen 
charge fluctuations and grain vibrations. We show that in the case of weak 
electromechanical coupling (i.e. when a change of the grain charge by e does 
not significantly affect the vibrational motion) the main results which can be 
obtained within the model introduced in Ref. [8] can be proven to be correct 
even for high resistance junctions. In this high ohmic limit a new scenario for 
the shuttle instability becomes possible leading to hysteresis in the current - 
voltage characteristics. 

The paper is organized as follows: In Section 2 we present our model composite 
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Fig. 1. (a) Simple model of a soft Coulomb blockade system in which a metallic 
grain (center) is linked to two electrodes by elastically deformable organic molecu- 
lar links, (b) Dynamic instabilities occur since in the presence of a sufficiently large 
bias voltage V the grain is accelerated by the same electrostatic force towards first 
one, then the other electrode. A cyclic change in direction is caused by the repeated 
"loading" of electrons near the negatively biased electrode and the subsequent "un- 
loading" of the same at the positively biased electrode. As a result the sign of the net 
grain charge alternates leading to an oscillatory grain motion and a novel "electron 
shuttle" mechanism for charge transport. 



Coulomb blockade system as well as a derivation of the equations governing 
the charge transport through the system. Then, in Section 3 we derive the 
dynamical equations for the case when the 'fast' grain oscillations have been 
averaged out so that only the 'slow' variation in oscillation amplitude remains. 
In Section 4 these equations are used to analyze the shuttle instability and 
the resulting current - voltage characteristics in our model Coulomb blockade 
system. The results are in good agreement with our earlier numerical results. 
The analytically solvable model is in addition useful for analyzing the nature 
of the loss of stability as the shuttle mechanism sets in. Such an analysis 
- supplemented by Monte Carlo simulations — is carried out in Section 5, 
where we find that depending on the resistances of the tunnel junctions the 
transition from the static regime to the shuttle regime can be associated either 
with smooth increase in the amplitude of the self-oscillations or with a jump 
in the amplitude at the transition. Adopting nomenclature from the theory 
of oscillations [9] we are dealing with either soft or hard excitation of self- 
oscillations, where hard excitation is associated with a hysteretic behavior of 
the current as the bias voltage is swept up and downQ. Finally, in Section 6 

1 In the language of phase transitions — taking the oscillation amplitude to be the 
order parameter — the 'soft' case corresponds to a second order transition and the 
'hard' case to a first order transition 



3 



we present our conclusions. 



2 Model System 

In this Section we present a model of the simplest possible composite Coulomb 
blockade system that retains the properties of interest for us. It consists of one 
small metallic grain of mass M connected by elastic molecular links to two 
bulk leads on either side, as shown in Fig. la. The electrostatic potential of 
the grain is a linear function of the grain charge Q and the bias voltage V 

where X denotes the displacement of the grain from the equilibrium position 
in the centre of the system. This relation defines the capacitance C(X) that 
appears in the theory. In a symmetric situation, as the one considered in this 
paper, 4>(Q, V, X) = 4>(Q, — V, —X) and hence the coefficients C(X) and a(X) 
are even and odd functions of X respectively. For small displacements we have 
4> « Q/C, where C = C(0). The magnitude of C corresponds typically to the 
size of the grain. There are three different forces acting on the grain; a linear 
elastic restoring force F e \ = —kX, a dissipative damping force Fd = —^dX 
and an electrostatic force F q = 

F= Ue -^(Q L -Q K ), 

where the electrostatic energy U e of the system and the charges in the left 
(right) lead Ql (Qr) are considered to be functions of Q and V. The function 
T is bilinear in Q and V hence 

For a symmetric junction 

F(Q,V,X)=F(Q,-V,-X) 

implying that $(X) is an odd function and that C(X) and b(X) are even. 
Then, for small displacements, 

F q = -Q, L- 1 ee M{X) 

L dX X=Q 
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and one can consider £ = V/L as an effective electrostatic field induced by 
the bias voltage acting on the grain between the leads. This way we arrive at 
the estimate L ~ distance between the leads. 

Knowing the forces that act on the grain, we may now introduce its equation 
of motion, 

MX + ld X + kX = £Q(t) . (1) 



If the Coulomb charging energy, U c = e 2 /C, satisfies U c 3> h/RC, f3~ l where R 
is the characteristic tunneling resistance and j3 is the inverse temperature, one 
can expect a strong quantization of the charge Q(t) in units of the elementary 
charge e, 

Q(t) = en(t) , 



where n(t) will be a step function which can take on only integer values. 
Changes in n(t) with time are due to quantum transitions of electrons between 
the grain and the leads. According to the 'orthodox' Coulomb blockade theory 
[2] the probability for the transition 

(n, Ql,r) -> (n ± 1, Ql,r =F e) 



to occur during a small time interval At when the grain is located at X can 
be expressed as 

W[%(n, X, At) = At — L^ rffi(n, X) , (2) 



where 



Here the function / is defined as 



±n--). (3) 



f{x) ~ 1 - eM-PU c x) ' (4) 



where i?L(R)(^) is the tunneling resistance of the left (right) junction which 
depends exponentially on the distance between the grain and the respective 
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lead. In this paper Rl(X) = Rr(— X) = Rexp(X/ A) where we refer to A 
as the tunneling length. Depending on the material of the reservoirs and the 
insulating links A can be estimated to lie within 0.05A— 3A for direct tunneling 
from the electrodes to the grain. Due to the strong exponential dependence 
on the tunneling resistances the variations in capacitance with position are 
relatively small and are therefore neglected. The equations (l)-(4) hence define 
our model system. 



3 Limit of Weak Electromechanical Coupling 



In this Section we describe how the model system introduced above can be 
solved analytically for the case of weak electromechanical coupling. In later 
Sections this analytical solution will be compared to 'exact' Monte-Carlo re- 
sults and found to be very useful in the further analysis of the model system. 
The key is to average over the fast grain oscillations so that we are left with a 
set of equations describing only the slow variations in the amplitude. We start 
by considering the typical scales of our parameters. It is natural to use A as a 
characteristic length scale. Furthermore the condition for observing Coulomb 
blockade effects requires that we operate with voltages V of the order of the 
Coulomb blockade offset voltage V = U c /e = e/C. Introducing the dimen- 
sionless variables x = X/X and v = V/Vq the equation of motion for the grain 
(1) turns into 

x + u 2 x = — vn{t) — 7X , (5) 



where uj = yk/M is the elastic oscillation frequency and 7 and f2 2 are defined 
through 7 = 7rf/M and 

02 = e2 = h.*ku? 

MLCX L E x ' 



Here E x = k\ 2 /2 is the energy of harmonic mechanical vibrations with ampli- 
tude A. For a nanoscale grain and typical organic junctions Q 2 /uj 2 = e ~ 10 -2 . 
Taking 7/cj to be of the same orderQwe can separate out slow variations in os- 
cillation amplitude by averaging over the fast oscillations. This is conveniently 



Although the analytical model presented in this paper is dependent on treating 
the electromechanical coupling and the damping as small perturbations to a simple 
harmonic oscillator, computer simulations reveal that the qualitative behavior of 
the system persists even when these terms are large. 
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done by looking for a solution to (5) of the form 

x(t) = x(t) sin(atf) , 



(6) 



where x(t) is a slowly varying function such that x(t) ~ euo. Substituting (6) 
into (5), multiplying by cos(u;t) and averaging over a time interval T yields 

dx Q 2 1 

dt = 2uj V ^ C0S ^W) ~ 2 7 ^ ' ^ 



The time average is indicated by brackets and defined as 

t+T/2 

(9(t)}, = ^ J 9(r)dr, 

t-T/2 



where the index x indicates that the averaging process is performed for har- 
monic oscillations with constant amplitude x. The value of T has to be chosen 
to obey the double inequality 

to' 1 , cur 1 = RC < T < t s = , 



where ts ~ 1CT 9 — 10 _8 s is the characteristic time for changes in the vibration 
amplitude x(t) and cu^ 1 ~ 10 _10 s is the typical time for charge redistribution 
between the leads and the grain. With this choice of T the grain will perform 
many oscillations which differ very little in amplitude during the averaging, 
i.e. 

5t(t + T)- 5t(t) 1 
x(t) 

Introducing the dimensionless mechanical vibration energy 

Muo 2 (\x(t)f 



E(t) 



2E. 



and making a change of variables in (7) we get 



H F O 2 

-p = v—W(E) — 7-E (8) 
dt uj 

W(E) = VE (n(t) cos(ujt)} E . (9) 
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The physical meaning of the term containing W(E) is that energy is pumped 
into the system when the grain is oscillating. Due to the correlation be- 
tween the charge fluctuations n(t) and the position of the grain, x(t) = 
\J E{t) sin(cjt), W(E) will be nonzero. In order to describe this correlation 
we introduce the correlation function 

P n (ip, E) =n (5 nMt) 5 (sin \{ut - </?)) ) e , (10) 



where the average is taken over a time T that allows the grain to perform 
N = Tuo/2-k complete harmonic oscillations with constant energy E. The 
grain will now pass through the particular point 

(E,ip) = (x,x) = (x E (ip),x E (ip)) = (VEsm(ip),ujVEcos(ip)j 



in a two-dimensional 'phase space' N times and P n (ip,E) will thus be the 
relative number of times the grain passes this point with charge Q = en. 
Using the definition (10) of P n ({p,E) we rewrite (9) as 



W(E) = ^Jd V cos(^)g(y2, E) (11) 
q( ( p,E) = '£nP n (<p,E). (12) 



In Appendix A it is shown how one can obtain a differential equation for 
P n ((p,E). The equations (8), (11) and (A. 6) then completely describe the 
behavior of the model and are stated below in their final form, 



^L = V ^L-W(E)- 1 E (13) 



W(E) = X_ J dy? cos(v?) £ nP n ( V , E) (14) 

dP [^ E) = u R Q(VE sin y?)i%, E) . (15) 
dip 

Here P((p, E) is a vector containing P n (<p, E) and the components of the matrix 
Q appearing in (15) are 

Qn, m (x E ((p)) = -6 nim [F%(n, (p) + T E (n, tp)\ + 5 n , m±1 T E {n, tp) 

r|(n, (p) = T^(n, x E (<p)) + Tg(n, x E (ip)) 
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and we have denned a dimensionless charge relaxation frequency, 

= • 

UO 

If we formally solve equations (13-15) any observable characterizing the system 
can be evaluated. One such observable is the average current through the left 
and rights leads respectively. In Appendix A this current is shown to be 

/l,r = ^ / dp£#*(>/Esinp)P B fo£), (16) 



where 



# Vising = u R e^ si ^ [f(v/2 T n- 1/2) - f(±n - v/2 - 1/2)] . 

(17) 

By using the definition (17) of the partial currents and Eq. (12) together with 
(A. 6) one comes to the expression 

= E sin ¥>) - ^(v^siny,)) P n (<p,E) . (18) 



This form of our equations is especially useful when v = (2n + 1) at zero 
temperature, i.e. when the function f(x) = x6{x) and when the voltage is 
chosen at the point where a new channel is about to switch on. The differential 
equation (18) for q then simplifies to 



dq(^E) 
dip 



= u R ( (1 — v) sinh(v / £' sirup) — 2q(ip, E) cosh(v / £' sirup) J . (19) 



Instead of n coupled differential equations for P n there is now only one for 
q(<p,E). 

We are now ready to apply our analytical solution of the model Coulomb 
blockade system to an analysis of the shuttle instability. This will be done in 
the next two Sections. 



4 Analytical and Numerical Analysis of the Shuttle Instability 



In this Section we will use both the approximate analytical approach of the 
previous Section and an 'exact' numerical scheme to analyze the shuttle in- 
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stability in our model Coulomb blockade system. Using first the analytical 
approach we conclude from Eq. (13) that the stationary regimes of our system 
are defined by 

n 2 

v—W{E) = >yE. 
to 



This equation always has one trivial solution E = 0. When this solution is 
stable it corresponds to the system being close to the point of mechanical 
equilibrium only subject to small deviations due to charge fluctuations. We 
will refer to this regime as the static regime. For our symmetric system the 
current-voltage characteristics in this regime show no pronounced Coulomb 
blockade structure even though there are peaks in the differential conductance 
due to the switching on of new channels at voltages v n = 2n + 1 [1]. 

In Appendix B it is shown that for small E we have 

W{E) = a(v)E + 0(E 2 ), a(v) > . 



Hence, according to (13) the point (x — 0, x — 0) will become unstable if 

dE 



dt 



tt 2 

> <^ v — a(v) > 7 , 



i.e. when more energy is pumped into the system than can be dissipated. We 
define the critical voltage v c as the voltage when 

v c —a{v c ) = 7 . 

This equation cannot be solved for v c in the general case, but if we know a(v) 
for some specific values V\ and v 2 of the bias voltage and if 

tt 2 tt 2 
vi — a(vi) < 7 < v 2 — a(v 2 ) , 



it follows that v c lies in the interval [t>i,t>2] since a(v) is a continuous function 
of v. For the special case v — v n = (2n + 1) at zero temperature we only 
need to consider the single equation (19). By performing successive partial 
integrations of (14) using (19) and approximating sinhx ~ x and cosh a; ~ 1 
one finds that 

_ v R (v n - 1) 
a[Vn) ~ 2(1 + 4*4) ' 
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Hence the critical voltage can be established within one Coulomb blockade 
voltage, 




70; (1 
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< 1 



(20) 



If this critical voltage is exceeded the oscillation amplitude will increase as the 
voltage is raised above v c . The answer to the natural question whether this 
increase of amplitude will saturate at some definite value or not is determined 
by the behavior of W(E) as E tends to infinity. To investigate the behavior 
of W(E) for large values of E it is convenient to consider the contributions to 
W(E) from three different regions in (x,x)- space; a "right", a "left" and a 
"central" region, see Fig. 2. The "right" region is defined as xe(<p) > xe(<£o) = 




Fig. 2. Division of the phase space into three regions: 'left', 'center' and 'right'. In 
the left and right regions the charge distribution on the grain is determined by the 
exchange with the corresponding lead. In the center region charge exchange with 
both leads contribute to the charge distribution. 

haE, where the rate of the charge exchange between the grain and the right 
lead dominates over the exchange with the left lead. Similarly the "left" region 
is defined as xe{^p) < ~xe{^o) where the converse is true. The third region 
is the central region where |x_b(</?)| < x E (<Po)- 111 a symmetric system the first 
two regions contribute equally to W(E) and we need only consider the "left" 
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region. For this case it is convenient to represent Eq. (15) for P(ip,E) in the 
form 

- c\P - 

g L P = -e(E)r R g( V )^--e 2 (E)g 2 ( V )g R P, (21) 

where 

Ql,r = e ±x Q(x)\ x>>1 , e(E) = exp(-x E ( Vo )) = E' 1 
g(<f) = exp(V / E'(sin(/9 + sirup)) < 1 , 

i-e. Gl,r describes the charge exchange with the corresponding lead when the 
grain is disconnected from the other lead. For large E when e < 1 one can 
develop a perturbation procedure to solve (21). For the 'left' region one finds 
the expansion 

P^)=P L -e 2 (E)g 2 g i : 1 g R P L + ..., 

where P L , which satisfies the equation g L P L = 0, describes the charge dis- 
tribution for a grain in thermal equilibrium with the left lead. Therefore the 
charge on the grain will saturate at the value q L = nP^ and as a consequence 
the leading term does not depend on the position of the grain. This is why 
the contribution from the left (and right) region to W(E) will decrase with 
increasing amplitude. The contribution from the central region is restricted 
by the upper limit 4xe(<Po) max{g(<^)}. All the above considerations yield the 
inequality 

W(E) < AlnE + B when E > 1 . 



This inequality implies that at large amplitude the rate of dissipation, ^E, 
always exceeds the rate at which energy is pumped into the system leading 
to a final stationary state with the grain oscillating with a finite amplitude 
E = E s h if v > v c . This oscillating regime we refer to as the shuttle regime. A 
more sophisticated treatment reveals that 

W(E) = 2q L \nE + 0(1) when E > 1 . 



This leads to the following estimation for E sh as a function of the bias voltage 
at v — v c at zero temperature when 2q L ~ v 

Q 2 v\ fn 2 v 2 \ 

E sh = In . 

^7 V ^7 
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To support the preceeding discussion we have performed numerical simulations 
of the system. Figure 3 shows the charge on the grain as a function of its 
position in the shuttle regime. The thin lines are the result of a Monte-Carlo 




Fig. 3. Comparison between 'exact' Monte Carlo simulations and the analytical 
approach. The graph shows the grain charge as a function of grain position in the 
shuttle regime at bias voltage V = lOVb (Vo is the theshold voltage for lifting the 
Coulomb blockade). The thin lines are from stochastic simulations and the thick 
smooth curve is the result of our analytical approach to the system. From this 
graph we see that the charging-decharging process takes place in the center of the 
system while the the charge on the grain saturates at it's maximum (minimum) 
value outside this region. 

simulation of the system when v — 10 and the thick lines show the results of 
our analytical approach. The graph reveals the true stochastic nature of the 
charge on the grain as well as the validity of our analytical solution. Moreover 
one sees that the charging and decharging of the system takes place in a limited 
region around the center of the system while the charge saturates as the grain 
approaches the respective lead. To understand this behavior one can integrate 
Eq. (18) to get the following relationship for the charge on the grain 

9(0) - g(TT) = / Y,3nPn{<p) dtp- f T,jnPn(v) d V 

i n i n 

(Recall that g(27r) = q(0) and that q(0) [q(tr)] is the charge on the grain as 
it passes the center position moving right [left]). Substituting this expression 
into (16) one can get the current in the form 
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(22) 



The last two terms determine the current between the grain and the most 
distant lead. We can think of this contribution as a tunnel current through 
the central cross section of the system. The first term corresponds to the 
current through this cross section when the grain passes the point x — 0. 
This current exists only because of the oscillatory motion of the grain. We 
refer to this mechanically mediated current as the shuttle current. At large 
amplitudes we can estimate the contribution from the last two terms in (22) 
to be of the order 1/y/E. At the same time it follows from a perturbative 



treatment analogous to the one carried out above to determine the behavior 
of W(E) for large E that 



At zero temperature we have also 



where the square brackets [..] denotes the integer part of the argument. Thus 
we find for this case that 



From (23) one finds that for large amplitude oscillations a pronounced step- 
like behavior in the current appears in the symmetric case in contrast to the 
ordinary static Coulomb blockade case. This behavior is also seen in numeri- 
cal simulations of the system. In Fig. 4 the current-voltage characteristics has 
been calculated. As the voltage is raised above v c the current departs from the 
current obtained in a static double junction and distinct steps appear. 



5 Soft and Hard Excitation of Self-Oscillations 

We have in the previous Sections shown that the stationary point (x — 0, x — 
0) will become unstable for certain voltages v > v c and that the system will 
reach a self-oscillating regime with well defined amplitude. In this Section 





(23) 
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2 4 6 8 10 

v 



Fig. 4. Current through the Coulomb blockade system in the static regime (V < V c ) 
and in the shuttle regime (V > V c ). As the system enters the shuttle regime the 
current (solid line) deviates from the current for a static double junction (dashed 
line). After the transition to the shuttle regime distinct steps as predicted by Eq. (23) 
can be seen in the current even though we are modelling a symmetric system. The 
current in this figure and subsequent ones is normalized to the frequency of harmonic 
oscillations to demonstrate that the step height in the shuttle regime is proportional 
to oj. In order to make a comparison with the current in the static regime this current 
has been scaled by the same factor. 

we now show that the transition from the static regime to the shuttle regime 
can be associated with either soft or hard excitation of self-oscillations. The 
terminology is adopted from the theory of oscillatons [9]); in the case of soft 
excitation of self-oscillations the amplitude increases smoothly from zero at 
the transition point, while the oscillation amplitude jumps to a finite value in 
the case of hard excitation of self-oscillations. The former case occurs when 
W"(0) < and the latter when W"(0) > 0. 

We start with the case when W"(0) < then for small E and v — v c <C v c we 
have 

W{E) = [a(v c ) + a'(v c )(v - v c )\E - (3E 2 a,(3>0. 

The energy at constant amplitude satisfies according to (13) 

v — [cx(v) — j3E] = 7 , 
ui 

where v c a(v c )Q 2 /u = 7. Hence, in the vicinity of the transition the amplitude 
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v^W(E) 




O 



E 



Fig. 5. Schematic energy diagram for the case W"(0) < 0. The graph shows W(E) 
for a fixed set of parameters for three different voltages along with the line jE. 
When v = v ' < v c all the energy pumped into the system can be dissipated and 
the stationary point O is stable. When v = v c this point becomes unstable and at 
voltage v = v' > v c the system will reach a limit cycle with oscillation amplitude 
oc \[E determined by the intersection point P. 

of oscillation will increase smoothly as 



The development of the instability can be understood from the diagrams in 
Fig. 5. The graph shows W(E) for a fixed set of parameters for three different 
voltages along with the line 'jE. When v < v c the dissipation is larger than 
the pumping of energy into the system and the stationary point O is stable. 
When v — v c this point becomes unstable and at voltage v > v c the system 
will reach a limit cycle with oscillation amplitude oc \fE determined by the 
intersection point P. 

When W"(0) > the shuttle instability develops in a completely different 
way. Consider the diagrams in Fig. 6; the graph shows W(E) for a fixed set 
of parameters for four different voltages along with the line 'jE. Consider now 
the system being located in O at a voltage v < v c \. In this case the system 
is in the static regime and exhibits the same behavior as an ordinary double 
junction. As the voltage is increased above v& a second stable stationary point 
P' appears but the system cannot reach this point since O is still stable. At 
v = v c2 , O becomes unstable and the system "jumps" from O to Q 2 . This 
instability we refer to as hard since the amplitude changes abruptly from E — 
to E — E2 as the voltage is raised above v C 2- Now consider the case when the 
system is originally in the stationary point P" and the voltage is lowered. At 
v < v c2 , O becomes stable but cannot be reached by the system until v has 




(24) 
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Fig. 6. Schematic energy diagram for the case W"(0) > 0. The graph shows W(E) 
for a fixed set of parameters for four different voltages along with the line *yE. 
When v < v c i only O will be a stable stationary point. At v = v c i a second unstable 
stationary point Qi appears. For v c \ < v = v' < v c i we have two coexisting stable 
points O and P' leading to the hysteretic behavior of the system discussed in the 
text. At v = v C 2 O will become unstable and the system is determined to be in the 
limit cycle with amplitude corresponding to energy E2 at the intersection Q2. As v 
is increased above v C 2, the only stable stationary point left is P" corresponding to 
a limit cycle with amplitude oc \[E. 

dropped to v — v cl . At v cl the point Qi becomes unstable and the system will 
"jump" to O. This transition is characterized by an abrupt drop in amplitude 
from Ei to E — at v — v&. Since v c \ < v C 2 the system will obviously exhibit 
a hysteretic behavior in the transition region. 

By using the simplified equation (19) for the charge valid at zero temperature 
and v — 2n + 1 one can develop a perturbation expansion for q((p) for small 
E. Solving the system to third order in \f~E~ we find an expression for W"(Q) 



From this result follows that the instability is associated with soft excitation of 
self-oscillations if tr = z/^ 1 < 2\f?> and with hard excitation of self-oscillations 
if r R = v K l > 2^3. 

Monte Carlo simulations of the system support the existence of two different 
types of instabilities. In Fig. 7 the result of a simulation for the case when 
tr < 2y / 3 is shown. As expected the amplitude grows as in (24). The transition 
to the shuttle regime is visible in the current - voltage characteristics as a 
lowering of the current compared to the current for a static double junction. 
This decrease can be understood from the analytical expression (23) for the 



'W(E) 
dE 2 



E=0 




(25) 
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Fig. 7. Soft excitation of self-oscillations from Monte Carlo simulations, (a) shows 
the root mean square amplitude as the voltage is raised above V c when tr < 2\f?>. 
The predicted square-root increase of amplitude (Cf. Eq. (24)) can clearly bee seen 
here. The slight increase in amplitude before the transition comes from fluctuations 
in the charge on the grain, (b) The loss of stability is reflected in the current - 
voltage characteristics. As the system enters the shuttle regime the current deviates 
from the current seen in ordinary static double junctions (dashed line). 

current . 

A similar simulation for a case when r R > 2\/3 is shown in Fig. 8. The pre- 




Fig. 8. Hard excitation of self-oscillations from Monte Carlo simulations, (a) shows 
the root mean square amplitude at different voltages. As the voltage is raised from 
V = there is a discrete jump in amplitude at V = V C 2- Lowering the voltage again 
from V = 3Vo to zero reveals the hysteretic behavior as the amplitude drops to 
zero at V = V&. (b) The hysteretic behaviour is clearly visible in the current. Since 
i^« 1 the current in the shuttle regime lies very close to the value / = euj/ir as 
predicted in (23). 

dieted hysteretic behaviour is clearly visible along with the expected quanti- 
zation of the current in the shuttle regime. 
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6 Conclusions 



We have analyzed by both numerical and analytical methods charge transport 
by a novel 'shuttle mechanism' through the model shown in Fig. 1 of a self- 
assembled composite Coulomb blockade system [8]. A dynamical instability 
was found to exist above a critical bias voltage V c , which depends on the 
junction resistances in the system. In this 'shuttle regime' there is a limit 
cycle in the position-charge plane for the grain (shuttle) motion as shown in 
Fig. 3. Above V c the current - voltage curve has a step-like structure, a type 
of Coulomb staircase, as shown in Fig. 4 even though we are modelling a 
symmetric double-junction system. 

The transition from the static regime, where the grain does not move, to the 
shuttle regime can either be associated with soft excitation of self-oscillations, 
i.e. a continuous square-root increase in oscillation amplitude above V c as 
shown in Fig. 7 or by hard excitation of self-oscillations (the terminology is 
from the theory of oscillations [9]) implying a sudden jump in oscillation am- 
plitude at an upper critical voltage V c2 . In the latter case the system behavior 
is hysteretic, as shown in Fig. 8; when the bias voltage is lowered the static 
regime is re-entered at a lower critical voltage V c \. 

Hard excitation of self-oscillation appears in our model for large junction resis- 
tances, tr = RC > 2\^3/uj. At the same time the numerical value of V c oc \fR 
for large resistance junctions and may therefore be quite a bit larger than 
the threshold voltage Vq for lifting the Coulomb blockade. The current in the 
high-voltage shuttle regime may furthermore be distinctly larger than in the 
static low-voltage regime. This is because it is governed mainly by the elas- 
tic vibration frequency of the grain and not by the rate of charge tunneling 
from the grain at rest in the center of the system, equally far from both leads. 
Hence in an experiment one might find a sudden increase of current at a quite 
large voltage V c , where the shuttle mechanism sets in, rather than at the lower 
threshold voltage V , where the Coulomb blockade is lifted. This observation 
may be relevant in connection with the recent experiment by Braun et al. [3] 
where hysteresis effects in the current - voltage curves and an anomalously 
large threshold voltage (if interpreted within conventional Coulomb blockade 
theory) were reported for a system showing very large resistance at low bias 
voltages. The results of this work show that it may be possible to explain 
this experiment within a Coulomb blockade theory which allows for a shut- 
tle instability of the type described here. In order to establish whether this 
mechanism is truly responsible for the behaviour of the system in [3] one 
has to study the experimental situation in great detail to properly estimate 
the mechanical parameters (elastic vibration frequency and damping) of the 
insulating material. 
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A Appendix A 



In this Appendix we derive the differential equation (15) for the quantity 
P n ((p, E) defined in (10) and the result (16) for the average current /L,R.through 
the left and right leads. 

Denote by N n (ip,E) = NP n (ip,E) the number of times the grain passes 
through the point (E, ip) carrying charge Q = en and let 
A r m „(<^, E, At) be the number of times the grain has charge Q = em as it 
passes through the point (E, <p + a; At) provided that it passed through (E, ip) 
at a time At earlier with charge Q = en, i.e. 

N mn (<p, E, At) = Ntt (5 m7n (t+At)S n , n (t)S (sin (±(u;i - <f)))) E (A- 1 ) 

It's now easy to see that 

N n { v + uAt, E) = J2 N nm {^ E, At) = ]T N nm { v , E, At) 

+ U n ( v ,E)- ]T N mn (<p,E,At)\ (A.2) 
Using the definition of P n (ip,E) in (A.2) yields 



u f . a j. 7ti\ \- N nm (ip,E, At) 

+ U(.*.E)-T N i^ t] P^ (A.3) 

For large N the ratio N nm (tp, E, At) / N m (tp, E) is the probability for the tran- 
sition (n, Ql,r) — ► (m, Ql,r T ( m ~ n )) during a small time interval At. But 
the probability for this kind of event was given in (2) 
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<x=±S=L,R ' u J 



+ E E Wl(n-al,x E (tp),^-)p n - al (tp,E) (A.4) 
We can write this as a first order differential equation if we introduce 

r l( n ; <P) = r t( n , xe(<p)) + rg(n, x £ (</?)) 

= e"^ si ^/(±| =F n - ±) + e^ sin ^/( T | =f n - ±) (A.5) 

and use the relations (2) and (3) 

^ = ~ (rj(n, *>) + ri(n, *>)) P n fo 

+ ^ (ri(n + 1, ¥»)P n+1 (^ £) + - 1, v?)P„-i(v?, £)) (A.6) 
The average current through the left (right) lead can formally be written as 

7 L,R = e (^ *E°WL*(n(t),X(t),At)j (A.7) 



Within our approximation we now consider this average when the grain is 
oscillating with a fixed amplitude \\^E and arrive at the equation 

/l,r = -^t /dv?E E ff ^R (^Av^sin^At) ± (<S n)n(t) <f(sin ±(wt - p)) 

^ J n a=± 

(A.8) 

By using the definition of the partial currents (17) one gets eq. (16). 



B Appendix B 



In this Appendix we prove that the function W(E) is a linear function of E 
for small E with a positive coefficient. 

For small E we can expand (15) around \fE = to obtain 

dP ^ E) = (So + V^sin^i + ...) P{<p,E) (B.l) 
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where 

dg(x) 



Go = G(0) and Gi - ^ 



(B.2) 

x=0 



Expanding P((p, E) in the same fashion 

P(<p, E) = + y/EP<U(<p) + ... (B.3) 

and inserting this Ansatz into (B.l) we find that P^ is the solution to the 
homogenous equation 

G P (0) = (B.4) 
The correction to first order in y/E, P^\(p) is the solution to the equation 

( TR ^ ~ i5(1)(v?) = ^ lPi0) SinV (B - 5) 
this is a standard differential equation which can be solved to yield 



P {1) (f) = - t , 2 , 2 Gi pi0) cos ^ + S f R 2 A 2 Gi p(0) sin ^ ( B - 6 ) 



I + vlGl I + v\G 







where / denotes the identity matrix and the matrix divisions are symbolic for 
the multiplication whith the inverse operator. Recalling the expression (14) 
for W(E) and introducing the vector n having components n we find when 
inserting our expansion that we have to first order in E 



W(E) = — [ d(pcos(ip)h- P(ip,E) 
2ir J 

= f (h ■ jJ^&pw) + Om = aE + Om (B.7) 

However, to get instability we require that a defined as 

a = i- j dy? cos(y?) £ nP^ (<p) (B.8) 

must obey a > 0. To show this consider the auxiliary equation 

(r R -^ - So) P (1 \r) = GiP {0) 9(r) (B.9) 
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where g(r) can be any non singular real valued function. We now define 

g (1) (r)=^nPW(r) (B.fO) 

n 

and the functional A{g(r)} as 

AMt)} = s~%*1 (B11) 

Note that we have as a special case 

a = ^A{g(r) = smr} (B.12) 

From the theory of linear response 

oo 

q^(r) = J X (r-r')g(r')dr' (B.13) 



where x( T ~ T ') is the response function. Using (B.13) in conjunction with 
(B.ll) one gets 

A{g{T)] ~ JMfW (B ' 14) 



where the Fourier transform is defined here as 

g{u) = J dre^(r) (B.15) 

Causality together with the requirement that q^ (r) be real give two conditions 
Im{xH} ^ for uj ^ and Im{x(-^)} = lm{- X (oo)} (B.16) 



Since x is °dd and nonzero at uj = we conclude that 

sgn ilm{uJx(co)}) = Const (B.17) 

However, the response function x(^) is independent of the particular form of 
g{r) and hence we have that 

sgn(a) = sgn(A{sinr}) = sgn (A{g(r)}) (B.18) 
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To see that a > consider the special choice of g(r) 



Irl > T 

9o(r) = { (B.19) 
-fi'osgn(r) \t\ < T 

where T uj/ur. This choice of g has a very simple physical meaning. At 
r < — T the grain is immobile in the centre of the system carrying charge 
q = 0. At r = —T it instantanously jumps to a location close to the right lead 
and waits there long enough to acquire a negative charge q = —go- It then 
moves instantanously to the left lead and there gets the charge q = +go and 
then moves back to the middle. This process yields 

A{g (r)} = 3g /2T > (B.20) 

Hence we have shown that for small E 

W(E) = aE + 0(E 2 ), a>0 (B.21) 
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